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Monte Carlo switching moves ( "perturbations" ) are defined between two or more clas- 
sical Hamiltonians sharing a common ground-state energy. The ratio of the density 
of states (DOS) of one system to that of another is related to the ensemble averages 
of the microcanonical acceptance probabilities of switching between these Hamiltoni- 
ans, analogously to the case of Bennett's acceptance ratio method for the canonical 
ensemble [C. H. Bennett, J. Comput. Phys., 22, 245 (1976)]. Thus, if the DOS of 
one of the systems is known, one obtains those of the others and, hence, the parti- 
tion functions. As a simple test case, the vapor pressure of an anharmonic Einstein 
crystal is computed, using the harmonic Einstein crystal as the reference system. For 
this case, compared to the adaptive Wang-Landau algorithm, the proposed method 
is considerably faster. As a further example of the algorithm, the energy dependence 
of the ratio of the DOS of the square-well and hard-sphere fluids is determined, and 
the temperature dependence of the heat capacity at constant volume of the system 
calculated and compared with canonical Metropolis Monte Carlo estimates. By its 
inherent perturbation approach, the algorithm is suitable for those particular cases 
where the properties of a related system are well known. 
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I. INTRODUCTION AND BACKGROUND 



A complete description of an isolated system at energy E is given by the phase space 
volume, 

Q{E) = iVpiV / m®(E - H(q,p)) (1) 

where q,p are 3iV-dimensional vectors stating the positions and momenta, respectively, of 
the N particles, Q(x) the Heaviside step function, h Planck's constant and H(q,p) the 
Hamiltonian. Through this quantity — or the closely related density of states (DOS) u(E) = 
dVt/dE — the connection with the entropy of classical thermodynamics is established as one 
of S oc lnfi(-E') (Hertz definition) or S oc \nu(E) (Planck definition). These two definitions 
are not mathematically equivalent and, as pointed out and discussed by Dunkel and Hilbert 
(see Ref . Q] and references cited therein) , there is disagreement in the literature as to which 
one is correct. However, these two definitions become numerically the same for large systems. 
Indeed, if the system is great enough in the number of its degrees of freedom, fluctuations 
in its kinetic energy will be vanishingly small, the potential energy distribution will be 
Boltzmannian and the system can be said to be at equilibrium at constant temperature, 
which is a desirable situation as it can be reproduced more readily in reality, for which the 
systems studied are generally large in this sense. Unfortunately, the size of systems that can 
be investigated by computer simulation may still be far from adequately approaching this 
limit. 

The most common solution to this problem of size is to couple the system, in the mathe- 
matical sense, to an infinitely large heat reservoir at constant temperature, thereby creating 
a formally infinite system. This system is governed by the canonical partition function 
(CPF), which can be expressed as 

/oo 
dEu(E)e~ E/kT (2) 
-oo 

where k is Boltzmann's constant, and T the absolute temperature. This, or a mathematically 
equivalent, route to the CPF has been exploited in numerical methods such as the reference 
system equilibration (RSE) method,^ the histogramming^! and multihistogramPS! methods, 
the histogram reweighting method,^ the Wang-Landau (WL) samplingP^ multicanonical 
methodspH^ transition matrix method^ 14 * 15 * or the nested sampling (NS) algorithmic^ 
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The CPF is directly related to the free energy by 

A(T, V) = -kT In Q(T, V) (3) 

and thus knowledge of it enables one to compute the temperature or volume, V, dependence 
of any desired thermodynamic property. Because the integrand co(E)e~ E ^ kT is a sharply 
peaked function in E, it is numerically an easier task to obtain the partition function at a 
specific temperature than to obtain the complete DOS. If one is interested in free energies 
only at one or a few specific temperatures, especially low ones, direct methods^^l to the 
free energy, i. e. the CPF at pre-defined temperature, will always be more efficient, simply 
because they have a smaller region of integration about which to worry. The DOS approach, 
on the other hand, is more powerful when a range of temperatures is of interest, especially 
in systems or models where the CPF has no volume dependence, e. g. lattice models. The 
width of the temperature interval of interest implicitly defines the width of corresponding 
energy interval [E_, E + ] that one needs to consider in a numerical search for u(E). However, 
for continuous potentials, u(E) approaches the known DOS of the ideal gas at high E, as 
the kinetic energy contributions will dominate the potential energy ones. In these cases, E + 
can be defined independently of any temperature. 

In the WL method, originally developed for discrete model systems but generalized to 
continuous Hamiltonians by later authors,^^ the DOS is computed through a random 
walk subject to importance sampling whose weights are iteratively adjusted in an attempt 
to make all energies equiprobable. The weights that achieve this are reciprocal to the DOS. 
The main drawback of the method is that the statistics to estimate the weights need to be 
gathered anew between each update of the weights, leading in some sense to a "duplication 
of efforts." In this connection, a notable effort to improve the convergence rate was recently 
published.^ 

In the NS method, this is circumvented by having the random walk subject to a weight 
function that acts on a non-uniform distribution of energy segments, concentrating on the low 
energy regions. The limits of each energy segment are calculated "on the fly," starting from 
high energies and proceeding downwards by cutting the lower segments in two. The limits 
are set from the condition that the probability of encountering a configuration belonging to 
a given segment be equal to a predetermined function of the depth of that segment in the 
partitioning tree. Once a limit is found, it stays fixed and is never subject to reevaluation. 
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This removes the "duplication of efforts" of the WL algorithm. The most problematic 
case for this algorithm is for potentials that exhibit large regions of infinite energy, as in 
for instance the square-well fluid, because then the energy partitioning scheme cannot be 
gradual. There is hence no benefit in using the simplification of the hard molecular core 
with this method. 

In the RSE method, on the other hand, the system is coupled to a finite heat bath with 
which it is allowed to exchange energy. The DOS of the heat bath is presumed known. The 
combined system is evolved according to the microcanonical probability distribution and the 
probability of the system of having different energies E < E tot is histogrammatically tracked 
and related, up to an i? tot -dependent factor, to the sought system DOS. The main drawback 
is that the factor can only be calculated precisely for very low energy and smooth potentials, 
and good statistics is only obtained in a narrow range of E below -Etot, meaning in practice 
that several simulations at different E tot have to be run. By careful considerations of the 
continuity of the DOS, the ^tot-dependent factor may be extrapolated to higher energies in 
the end. 

In this Paper, I propose an alternative method: a new route to obtaining oj(E) assuming, 
like the RSE method, that the DOS of a different system is known. Unlike the RSE method, 
however, the idea is that the other system is also similar, and thus knowledge of its DOS 
is able to speed up the calculation by focusing on the difference between the two systems. 
This lessens the need for importance sampling. We shall develop the method in the next 
Section, and after that I will present some simple numerical examples of its use. 

II. DESCRIPTION OF THE ALGORITHM 

Consider two systems, for simplicity labeled as and 1, for which the DOS are ujq(E) 
(known) and u>\(E) (sought). Classical microcanonical sampling of either of these systems 
can be carried out efficiently if the potential energy depends on the configuration only. Under 
this restriction, one simulates a Markov process in configuration space using the acceptance 



if E > Uf and zero otherwise, where U{ and Uf denote the potential energies before and 
after, respectively, an unbiased trial move in configuration space. This equation represents 
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probability,™^ 





the ratio between the densities of the kinetic energy states for a M- dimensional configuration 
space with iV molecules, and is the proper weight function to use for the microcanonical 
ensemble where all accessible states are considered equiprobable. 

Let us now suppose that systems and 1 are "similar" in the sense that they share the 
same configuration space. Then system 1, differing only by its potential energy expression, 
can be viewed as the result of a perturbation on system 0. For instance, let us define the 
generalized Hamiltonian, 

H x (q,p) = H (q,p) + \U'(q) (5) 

where U'(q) is the perturbation depending on q only, and A an interpolating factor between 
the reference (A = 0) and fully perturbed (A = 1) system. Correspondingly, we may define 
£l\(E) according to eq. ([!]) after inserting eq. ([5]). Let us now consider the superensemble 
that includes A as an extra coordinate. We write its phase space volume as, 

Q(E) oc J dgdpdAdCe(E - H Q (q,p) - XU'(q) - ( 2 /2r]) (6) 

where rj is a formal mass associated with the A-motion and ( is a dummy variable of integra- 
tion for the formal momentum of this motion. The total energy of this ensemble is E which 
is different from the regular total energy E because it also contains a kinetic contribution 
r]X 2 /2 in addition to that of the regular coordinates. We are not interested in the ensemble 
at energy E. Therefore — without loss of generality — we let rj — >■ whereby E — >■ E and then 
take the derivative with respect to E and thus have 

Q(E) oc J dXio x (E) (7) 

Because this integral has a weighting factor which is A-independent, the acceptance proba- 
bility of each state is still only proportional to the density of its kinetic energy states, and we 
see that the random walk also along the A-coordinate should be governed by an unmodified 
eq. Q. 

We now define an "equilibrium constant" Kij(E) as the ratio between the number of 
cycles the Markov chain sampling the superensemble visits system A = Aj to the number of 
cycles it visits system A = Aj. The existence of this equilibrium constant is guaranteed by the 
detailed balance condition that the Markov chain fulfills. From the direct proportionality 
between the microcanonical probability and the DOS it follows that 

K rm- - <p « )i m 

K -' {E] -^m-m, (8) 
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where Pij is the acceptance probability for changing from A = Aj to A = \j, and (...)* denotes 
a microcanonical ensemble average over the system A = Aj. The last equality follows from 
the flux balance at equilibrium: to wit that the gross flux between two states, given by 
the product of the acceptance probability and the occupation probability, is equal but of 
opposing sign in the reverse directions, i. e. 

P E {U U U f )(E - uy 1 ^ 2 - 1 = P E (U f , U % ){E - Uf)^ 2 - 1 (9) 

where Ui and Uf denote the potential energies of two states defined by any two arbitrary 
sets of the non-momentum coordinates of the Hamiltonian. The quantity (E — uY' IN l 2 ~ l 
is proportional to the kinetic energy DOS, and because the configuration space is sampled 
subject to this bias, it is thus also proportional to the probability of being in a state of 
potential energy U. In integrating the equation over all configuration space, we identify the 
average of Pe on both sides weighted by this probability (up to a constant of proportionality), 
whence eq. (|8]) is obtained. It is at this point appropriate to stress that in the case when 
MN/2 = 1, eq. ^ does not hold and the algorithm, as here outlined, is not applicable. 
Such is the case of a single particle (N = 1) confined to two spatial dimensions (M = 2); it 
is never the case in three spatial dimensions (M = 3). Excepting that special case, we have 
that 

K m (E)= i U 1 K l , i+1 = ( ^ = ^§ i (10) 

and it is from this relation that oji(E) may be extracted, if ojq(E) is known, in addition to 
the ratio u\(E) /ojq(E) which is always obtained. The method outlined may be regarded 
as a special case of Bennett's method,^ but applied to the microcanonical ensemble. Al- 
ternatively, if one does not sample the transition probabilities, but instead propagates the 
system between the two states, it may also be regarded of the expanded ensemble 

method^ applied to the microcanonical ensemble; this, however, is a line of attack which 
we shall not pursue. 

Implicit in the derivation so far is that the minimum value of the potential energy ex- 
pression is to be independent of A. In other words, the "potential energy" of a configuration 
(q, A) is to be understood as the potential energy difference with respect to the global po- 
tential minimum over q keeping A constant. This follows from eq. ^ (unless we make Pe 
explicitly A-dependent) and the following argument. In eq. ([9]), the quantity E — U is the 
kinetic energy. Let us say that the maximum kinetic energy is E, obtained when U = Uq, 



which is hence the minimum potential energy. As the A-coordinate does not affect the kinetic 
energy, the potential energy U = U should correspond to the maximum kinetic energy E 
regardless of the value of A. This does not restrict the method in any formal sense but it may 
pose a practical hurdle for very complicated Hamiltonians for which energy minimization is 
difficult. This is especially true if several intermediate A- values are considered over a chain 
of gradual perturbations, if these affect the energy minimum in a non-trivial way. 

III. NUMERICAL EXAMPLES 

In principle, any two Hamiltonians for which we can calculate the requisite ensemble aver- 
ages (-Poi)o an d (-Pio)i can be used in this method. The algorithm hence does not pose any 
greater programming challenges than that of regular Metropolis Monte Carlo techniques. 
Here we shall consider two cases: the square-well fluid and a class of anharmonic Einstein 
crystals. The simplicity of the latter model is motivated by a desire to keep the computa- 
tional demands low as comparison with the computationally demanding WL algorithm is 
prohibitively expansive otherwise; the generalization to more degrees of freedom is trivial in 
all other respects. The square-well fluid, on the other hand, presents an interesting test case 
in that appreciable regions of its configuration space are of infinite potential energy. Last 
but not least, the simplicity of the numerical examples makes it easier to gauge the correct 
behavior to be exhibited by the calculation. 

Except where noted otherwise, all of the numerical examples to be presented have been 
compiled using the GNU C Compiler (version 4.4.3) and executed on a single 2 GHz core of 
the author's "Intel Core 2 Duo" laptop computer. Memory demands of the calculations are 
all insignificant. 

A. Anharmonic Einstein crystal 

I present in this section some simple numerical test cases on a class of anharmonic Ein- 
stein (AE) crystals. By this I mean crystals for which the CPF of N molecules can be 
written Q = q(T) 3N , where q(T) is the partition function of a single one-dimensional oscil- 
lator. The natural reference system to use when approaching the AE crystal is that of the 
harmonic Einstein (HE) crystal since q(T) for an harmonic oscillator is known analytically: 
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its (classical) DOS, 

is independent of energy. In this equation, m is the mass and kf the force constant. The 
anharmonic systems studied were governed by the potential energy expression 

u x (x) = x 2 + Xx 4 (12) 

where the term Xx A is the perturbation and x the displacement from equilibrium. The 
energy range was discretized into intervals of O.lwo(l), an d the calculation at each energy 
level proceeded until the relative change in u(E) between intervals of 9 x 10 7 cycles was less 
than e = exp(lCT 5 ) — 1. With this convergence criterion, the number of cycles that each 
energy level was sampled varied near-stochastically between 60 and 160 million with no clear 
trend for either A or E. The system was switched between the reference and perturbed states 
at intervals of 4.5 x 10 7 cycles. Each cycle consisted of either (80% probability) performing 
a microcanonical Monte Carlo displacement or (20% probability) sampling averages. 

Visual inspection (Figure [T]) reveals that the logarithms of the resulting DOS are quasi- 
linear in energy with a constant of proportionality directly proportional to A. The vapor is 
assumed ideal, so that the vapor pressure is given by 

kT 

Pvap(T) = (Xqf (13) 
where A is the thermal de Broglie wavelength. This equation is derived from the equality 
between the chemical potentials of the gas and crystal phases in the limit of iV — > oo. The 
results of these calculations are shown in Figure [2] plotted against temperature. They are 
useful as a "yardstick" of how large the perturbations considered here are in relation to real 
systems. 

We now turn to a comparison with the WL algorithm, keeping in mind that any compar- 
ison in terms of efficacy is not straightforward because the performance of WL algorithm 
depends non-trivially on how frequently and accurately one updates the WL weighting func- 
tion in the Markov chain, as well as on parameters like step size and so on. Still, even when 
restricted to just 17 energy subdivisions, I find that the perturbation approach runs at least 
15 times faster for these systems, and a smooth curve is difficult to obtain. In the numerical 
implementation, the energy histogram collected for each iteration of the WL algorithm was 
considered "flat" when the largest and smallest values differed by less than 10%, tested at 
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intervals of 10 6 cycles. In this example, the WL algorithm was initialized from the uniform 
DOS of the harmonic oscillator and so benefits in equal measure as the perturbation method 
from the similarity between the reference and perturbed systems.^ Results from the WL 
simulations are shown with opaque symbols in Figure [TJ As is evident — despite on average 
longer runs — the statistical noise is very much greater. 



B. Square- well fluid tetradecamer 

In this numerical example, we consider the square-well fluid as a perturbation of its hard- 
sphere analog. The DOS of the hard-sphere gas obeys the form u(E) = £(N, V)E 3N / 2 ~ 1 , 
where £(N, V) is an unknown function of the number of particles and volume. Our ignorance 
of the precise form of this constant of proportionality means that only the ratio between 
the DOS will be possible to provide completely. The unit of energy and temperature that 
we will use for the remainder of this section is the magnitude of the pair potential at unit 
distance and A = 1. The unit of length is the hard-core diameter. 

To be precise, the reference system interacts through the pair potential 




Uo(r) = { (14) 



where r denotes the intermolecular separation. The perturbation we introduce is 

!-l r < a 
(15) 
r > a 

so that the total pair interaction is written 

u\(r) = u (r) + Xu(r) (16) 

It is easy to realize, that for our chosen combination of systems, (-Pio)i = 1 and so there 
is, unlike in the previous section, no need to consider two ensembles explicitly. Thus, only 
configurations of the reference hard-sphere system have to be generated, and furthermore, 
these configurations are independent of E so that all averages {Pq\(E))q can be sampled 
simultaneously for a given density. In the calculations to follow, a random 95% of the Monte 
Carlo cycles consisted of propagating the microcanonical Markov chain and the remaining 
5% of accumulating averages. 



We let a = a/2, an arbitrary choice based purely on aesthethic appeal: it is the lattice 
constant of the close-packed cubic crystal. For this system, the energy minimum is — 52A 
for 14 molecules. Thus, we consider the total energy expression, 

U x ({ rij }) = 52A + £ u x {n 3 ) (17) 

i>j=l 

whose zero-level is independent of A. In the preceding equation, {r^} is the ordered set of 
all pairwise distances between the fourteen molecules. To satisfy the requirements of the 
microcanonical ensemble, we introduce the constraint that the cluster is confined to a fixed 
spherical volume, arbitrarily chosen to be either 5007r/3, i. e. corresponding to a radius of 
5, and a volume fraction of 1.4% ("low density"); or 1087r/3, corresponding to a radius of 3, 
and a volume fraction of 7/108 ("low-medium density"); or 97r/2, corresponding to a radius 
of 3/2 and a volume fraction of 14/27 ("high density"). 

For the propagation of the hard-core Markov chain, one molecule was moved at a time. 
At the volume fraction of 1.4%, the maximum step length was 3.0; at the volume fraction 
of 7/108, it was 1.0; and at the volume fraction of 14/27, it was 0.15. These maximum step 
lengths led to acceptance rates of 51%, 54% and 52%, respectively. The DOS was sampled 
in energy intervals of 1.4, starting at E = 32 for the low-medium density and covering 60 
different energy grid points. The calculations proceeded for at least 2 x 10 8 cycles, which on 
the author's machine took a little less than 3 minutes of real time, but considerably longer 
runs were found necessary to achieve good convergence in the lowest energy regions, where 
up to 20 minutes could be necessary. A refined attack would distribute the energy grid 
unequally over the energy range. 

One interesting aspect of the way we have defined the perturbation in A is the approximate 
self-similarity that arises because of its linearity. For large values of E/X, we have that (U\) oc 
A, regardless of which ensemble the average is taken over, meaning that uj\/ a (E) = U\(aE), 
where a is any positive real number. This breaks down when E/X <C 1, because the linearity 
does not hold anymore. Nevertheless, by letting a — > oo we see through this formula that 
uj\{E) — » uq(E), when E — > oo. An indication that the computer code is well and working is 
that the DOS for the square-well tetradecamer (given on the logarithmic scale with respect 
to the reference system in Figure |3| actually shows this mathematically proven convergence 
on that of the hard-sphere tetradecamer at high energies. The algorithm also runs quicker 
until convergence in those cases. The interesting part, where convergence is also a bit more 
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problematic, is for the low energy regions where the DOS of the square-well tetradecamer 
exhibits a clear deviation from its hard-sphere counterpart. The depth of this "dip" in the 
curve is decreased when the density is increased. It is easy to see why this should be by 
considering the close-packed density where the molecules have no liberty of movement left, 
conditions under which the hard-sphere and the square-well fluid are indistinguishable. 

In Figure |4j I show the temperature dependence of the constant- volume heat capacity of 
the coupled system at the "low-medium" density corresponding to the volume fraction of 
7/108. The heat capacity was calculated through the statistical mechanical relation 

, d\nQ , 9 <9 2 lnQ 
c v = 2kT ^r + kT ^fr- ( 18 ) 

The broad peak in this function at around T « 0.9 is characteristic of a first-order phase 
transition far from the thermodynamic limit in contradistinction with the singularity that 
one obtains for the infinite system, even if contrary to the case of Ref. [36] it is clear from the 
density in this case that it is question of a gas-liquid rather than a liquid-solid transition. 
At high temperatures, we expect the translational equipartition value of C v = 3N/2 = 21 
to hold, and this is borne out by the graph. Moreover, this is also the limiting heat capacity 
at low temperatures, since the law of Dulong and Petit does not hold at for the square- 
well fluid. This is because the potential is not analytical, and so there is no first-order 
quadratic potential energy term to contribute to the heat capacity. This gives rise to a 
largely symmetric peak in the heat capacity. For comparison, the heat capacity calculated 
from regular constant- volume Monte Carlo simulations and the fluctuation formula, 

are also shown in Figure |4j It is to be noted that these simulations are very difficult to 
converge in the low-temperature regime, not the least because of the numerical instability 
that arises from the T 2 denominator for small T. 



IV. CONCLUDING DISCUSSION 



In this Paper, it has been shown that calculating the CPF through the DOS by a pertur- 
bation method is a viable alternative to other techniques if the DOS of a related system is 
known. Compared to the exact WL method, the present alternative is considerably faster, 
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at least for the simple case considered. Technically, the algorithm amounts to sampling (at 
most) two microcanonical ensemble averages and so must be considered very simple. Indeed, 
one would only need to add a couple of lines of code to pre-existing molecular dynamics pro- 
grams, for instance, to implement this algorithm; and it would require also but very modest 
modifications to most Monte Carlo programs to implement the microcanonical average. The 
greatest obstacle to a painfree implementation of this method is that the potential energy 
minimum value has to be independent of A, requiring at the very least that efficient en- 
ergy minimization can be carried out on the systems of interest. A mitigating factor is 
the obvious fact that for any method or algorithm to calculate the low-energy DOS, such 
energy minimization must be carried out implicitly. Systems for which energy minimization 
is difficult, for whatever reason, are thus inherently difficult cases to calculate the complete 
DOS by any method. Incidentally, efficient energy minimization is also a prerequisite of the 
WL-like algorithm of Soudan et a/P^ 

One way to remove the need for a prior energy minimization is to employ a path between 
the systems of interest which is independent of the potential energy minimum. The most 
direct way of achieving this is — in the absence of periodic boundary conditions — a scaling of 
the volume or, equivalently, the intermolecular lengths between different densities. Moreover, 
in the low density limit the DOS is known as that of the ideal gas, so there is a natural 
reference system to use. An added benefit of scaling over the length scale, instead of the 
energy scale, is that the DOS is obtained as a function of both the volume and energy. 
Nevertheless, such calculations would come at considerable extra cost because the steep 
intermolecular repulsion means that very gradual perturbations have to be considered. 

One advantage of the perturbation method is its ability to calculate u(E) at any E- value, 
independently of the E-range one ultimately considers, which means that the DOS can be 
gradually accrued from completely separate simulations without the need of having to de- 
cide on a discretization scheme beforehand. This also opens up a vast array of possibilities 
for further improvement. For instance, a "smart", e. g. automatic and non-uniform, dis- 
cretization of the energy levels when calculating the DOS, so that those regions where the 
DOS varies most rapidly are sampled most thoroughly, is a natural extension, somewhat 
analogous to the energy segment partitioning of the NS method. 

However, the greatest drawback of the method is that prior knowledge of the DOS is 
generally very scarce. This limits the optimal applicability of this method because the 
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repertoire of systems with known DOS does not necessarily include those that are related 
to the system of study. It is therefore foreseeable that this algorithm will be most useful 
in conjunction with another method to calculate the DOS. Like this, once obtained for one 
system, a whole series of related systems will be amenable to structured investigation. Also, 
unless the absolute DOS is needed (to compute, for instance, a phase diagram) in many 
situations entropic differences may suffice. 
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FIG. 1. Logarithm of the DOS of oscillators subject to the potential energy of eq. (12) for A = 
0.1,0.2,0.5 given with respect to the A = 0.0 reference DOS. Transparent symbols are calculated 
according to the perturbation method. Opaque symbols are from 17 iterations of WL sampling 
over 17 energy intervals. The points at E = 0.1 are made to coincide by design as the WL sampling 
does only provide the DOS up to an undetermined factor that is unique to each run. 
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FIG. 2. Deviation of the vapor pressure of the AE crystal from the HE crystal reference. 
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FIG. 3. Difference between the logarithms of the DOS of the square-well and hard-sphere tetrade- 
camers. 
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FIG. 4. Constant-volume heat capacity of the square- well tetradecamer as a function of temper- 
ature at a volume fraction of 7 / 108. The broad peak in the heat capacity is indicative of a 
first-order phase transition of a finite-sized system. The dashed line is the translational equiparti- 
tion heat capacity of 21, as well as the heat capacity of the hard-sphere reference system. Squares 



denote heat capacities calculated from canonical Monte Carlo simulations according to eq. (19). 
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